最大概似估計法介紹(1)

• 在我們學過生物統計學後，我們學過t test的求解公式，也學過linear regression的求解公式，但你是否知道Logistic regression的求解公式呢?

– 讓我們先複習一下linear regression的求解公式推導過程：

1. 先定義我們的預測式(y = b0 + b1 * x)

2. 再定義我們的求解目標(希望讓殘差平方最小化)

3. 偏微分求解殘差平方方程(了解極小值出現的位置)

• 我們先省略第三個部分(通常你們會覺得第三個部分最難)，你們是否注意到前兩個部分是自行定義的?

最大概似估計法介紹(2)

• 讓我們先忘掉線性迴歸及最小平方法(最小平方法是18世紀發展出來的)，在這裡我們先回到統計學的原始公理(統計學是19世紀末的學科)，想想統計學如何描述母群以及樣本。

• 我們先用一個最簡單的例子，假設母群疾病盛行率為p，而我們對母群做一次抽樣(樣本數 = 80)，我們在這次抽樣中抽出了49個有病的人，31個沒有病的人…

– 先考慮一個比較簡單的問題，假定台北該疾病盛行率為1/3，台中盛行率為1/2，高雄盛行率為2/3，請問我們這次的抽樣最有可能是在哪個城市呢?

• 統計學在描述母群以及樣本時，是在考慮『試問母群參數為何(母群參數未知)，則抽出此樣本之機率最大』，我們的任何推論都以此為核心。

– 當我們思考到這步時，試著算出在台北抽出49個有病及31個沒病的機率吧(請用二項分布計算，但注意二項分布的基本假設)：

Taipei.p = choose(80, 49) * (1/3)^49 * (1 - 1/3)^31
print(Taipei.p)

## [1] 2.078886e-07

– 接著算出在台中及高雄抽樣的樣本機率吧吧：

Taichung.p = choose(80, 49) * (1/2)^49 * (1 - 1/2)^31
print(Taichung.p)

## [1] 0.0118359

Kaohsiung.p = choose(80, 49) * (2/3)^49 * (1 - 2/3)^31
print(Kaohsiung.p)

## [1] 0.05449674

• 請問我們的樣本最有可能在哪個城市所抽出的?

最大概似估計法介紹(3)

• 同樣的問題，假設我們現在很確定樣本是從台北抽出的，但我們事先不知道台北的疾病盛行率，故我們現在不知道p為多少，現在該如何求解呢?

sample.p = function (p) {
  choose(80, 49) * (p)^49 * (1 - p)^31
}
sample.p(0.1)

## [1] 5.459073e-29

sample.p(0.3)

## [1] 5.402339e-09

sample.p(0.5)

## [1] 0.0118359

sample.p(0.7)

## [1] 0.02270722

sample.p(0.9)

## [1] 8.193775e-12

• 核心概念就是，找到一個p使抽出此樣本的機率最大化，而這個概念就是『最大概似估計法』

– 有了這樣的概念後，該怎樣求解呢?

最大概似估計法介紹(4)

• 以剛剛的問題為例，求解的方法有非常多種，最直覺的方式是把0到1中間的數字都帶公式算一遍，如我們利用迴圈來求解(精確度至小數第三位)：

seq.p = seq(0, 1, by = 0.001)
result = numeric(length(seq.p))
for (i in 1:length(seq.p)) {
  result[i] = sample.p(seq.p[i])
}

which.max(result)

## [1] 613

seq.p[which.max(result)]

## [1] 0.612

• 看來結果是0.612！我們把圖畫出來看看：

plot(seq.p, result, type = "l", xlab = "母群參數", ylab = "樣本機率")

最大概似估計法介紹(5)

• 當然，數學家不會用這麼沒有效率的方式求解，數學家在這個問題上會定義公式並且求解：

• 結果發現，當p為0, 1, 49/80時有極值，但0, 1是出現極小值，而49/80是極大值，故答案為49/80

• 我們知道大家數學都不好，所以要解這種問題我們在R裡面可以請他幫我們解。

– 我們需要使用套件「stats4」內的函數「mle」(但她只能求最小值出現的位置，不能求最大值，所以要改寫我們的sample.p函數)：

– 請注意，函數「mle」並非使用數學上的微分求解，他使用的方式我們會在下一節課再詳細介紹！

library(stats4)
sample.p = function (p) {
  -choose(80, 49) * (p)^49 * (1 - p)^31
}
fit = mle(sample.p, start = list(p = 0.5), method = "SANN")
fit

## 
## Call:
## mle(minuslogl = sample.p, start = list(p = 0.5), method = "SANN")
## 
## Coefficients:
##         p 
## 0.6124266

49/80

## [1] 0.6125

• 答案準嗎?

練習-1

• 這邊有一串數列，假定他們是從一個常態分佈的母群抽出來的，我們現在不知道她的mean以及sd，請嘗試用最大概似估計法求解：

x = c(1, 7, 5, 6, 8, 3, 2, 9, 4, 5, 3)

– 你可能不清楚該怎麼求得常態分佈的機率，可以試著使用函數「dnorm」：

dnorm(0, mean = 0, sd = 1)

## [1] 0.3989423

dnorm(0, mean = 0, sd = 2)

## [1] 0.1994711

dnorm(0, mean = 1, sd = 2)

## [1] 0.1760327

dnorm(1, mean = 0, sd = 2)

## [1] 0.1760327

– 請與直接計算做比較

mean(x)

## [1] 4.818182

sd(x)

## [1] 2.522625

• 你發現了什麼?

以最大概似估計法求解線性迴歸(1)

• 現在我們要用最大概似估計法求解線性迴歸了，但在開始前必須下兩個定義：

1. 先定義我們的預測式(y = b0 + b1 * x)

2. 再定義我們的求解目標(希望讓樣本機率最大化)

• 樣本機率怎麼求呢?我們可以「假設」殘差為某個不特定的常態分布(但平均數必須為0)，接著我們就能帶入任意b0及b1求得樣本機率了！

– 這裡要注意一點，因為抽到每個個案的機率相當低，之後還要累乘會變得更小，可能會小到電腦無法紀錄，因此把它做對數轉換能有效解決這個問題！

– 另外，本來機率是累乘，現在取完對數後變成累加！

x = c(1, 2, 3, 4, 5)
y = c(6, 7, 9, 8, 10)

linear.p = function(b0, b1) {
  y.hat = b0 + b1 * x
  res = y - y.hat
  mean.res = 0
  sd.res = sd(res)
  log_p.res = dnorm(res, mean = mean.res, sd = sd.res, log = TRUE)
  return(-sum(log_p.res))
}

fit1 = mle(linear.p, start = list(b0 = 0, b1 = 0), method = "SANN")
fit1

## 
## Call:
## mle(minuslogl = linear.p, start = list(b0 = 0, b1 = 0), method = "SANN")
## 
## Coefficients:
##        b0        b1 
## 5.3000118 0.8987946

以最大概似估計法求解線性迴歸(2)

• 讓我們把結果跟傳統代公式算的結果做比較吧！

fit1

## 
## Call:
## mle(minuslogl = linear.p, start = list(b0 = 0, b1 = 0), method = "SANN")
## 
## Coefficients:
##        b0        b1 
## 5.3000118 0.8987946

fit2 = lm(y~x)
fit2

## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##         5.3          0.9

• 結果是不是很像?

以最大概似估計法求解線性迴歸(3)

• 統計分析中除了係數外，標準誤也相當重要(這樣才能求得p value)，因此我們也要比較一下他們的變異數-共變異數矩陣：

vcov(fit2)

##             (Intercept)           x
## (Intercept)   0.6966667 -0.19000000
## x            -0.1900000  0.06333333

• 這裡要注意，由於fit1是「S4 class」，這裡是我們課程中第一次遇到，注意拆解它的手段：

slotNames(fit1)

## [1] "call"      "coef"      "fullcoef"  "vcov"      "min"       "details"  
## [7] "minuslogl" "nobs"      "method"

slot(fit1, "vcov")

##            b0          b1
## b0  0.4370205 -0.11400892
## b1 -0.1140089  0.03800375

fit1@vcov

##            b0          b1
## b0  0.4370205 -0.11400892
## b1 -0.1140089  0.03800375

• 好像有點落差！

以最大概似估計法求解線性迴歸(4)

• 這其實是因為標準誤的底數是N或是N-2的關係，樣本夠大就不會有影響了

x = rnorm(1000)
y = 3 + 2 * x + rnorm(1000)

fit1 = mle(linear.p, start = list(b0 = 0, b1 = 0), method = "SANN")
fit2 = lm(y~x)

fit1@coef

##       b0       b1 
## 3.034672 1.981693

fit2$coefficients

## (Intercept)           x 
##    3.036488    1.983512

fit1@vcov

##               b0            b1
## b0  9.528576e-04 -7.396710e-06
## b1 -7.396710e-06  1.004551e-03

vcov(fit2)

##               (Intercept)             x
## (Intercept)  9.538094e-04 -7.418096e-06
## x           -7.418096e-06  1.006548e-03

練習-2

• 現在輪到你實作邏輯斯迴歸了，先讓我們看看他們的關係式：

– 可以轉換成這樣：

• 所以我們該怎樣計算每個樣本的機率呢?

x = 0:10
y = c(0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1)

b0 = -3
b1 = 0.5
p = exp(b0 + b1 * x)/(1 + exp(b0 + b1 * x))
p

##  [1] 0.04742587 0.07585818 0.11920292 0.18242552 0.26894142 0.37754067
##  [7] 0.50000000 0.62245933 0.73105858 0.81757448 0.88079708

• 接著要怎樣算整體機率呢?舉例來說以第一個數值為例，抽中的機率只有0.0474，因此若y為1則機率了0.0474，若y為0則機率為1-0.0474，之後請你寫出累積機率函數，並使用函數「mle」求解！

– 請與函數「glm」比較結果：

fit3 = glm(y~x, family = "binomial")
fit3

## 
## Call:  glm(formula = y ~ x, family = "binomial")
## 
## Coefficients:
## (Intercept)            x  
##     -1.4719       0.3409  
## 
## Degrees of Freedom: 10 Total (i.e. Null);  9 Residual
## Null Deviance:       15.16 
## Residual Deviance: 12.64     AIC: 16.64

自訂關係式(1)

• 我們現在很清楚若我們想要使用最大概似估計法的條件，在開始前必須下兩個定義：

1. 先定義我們的預測式

2. 再定義我們的求解目標

• 讓我們來嘗試自己製造一個方程式，並試著求解：

– 請到這裡下載範例資料。

dat = read.csv("data/Patient2.csv")
head(dat)

##            t        Y
## 1 0.00000000 232.9187
## 2 0.08333333 234.6244
## 3 0.16666667 234.1379
## 4 0.25000000 233.6271
## 5 0.33333333 232.8027
## 6 0.41666667 231.9975

• 這筆資料是一個經過放射碘治療的病患，在病房裡隨著時間身體放射性物質殘留的變化，其中t是時間，Y是感應器所偵測到的值

– 根據我們的醫學及物理學知識，我們了解到時間與放射物質殘留量的關係應如下所示：

• 其中，

D_1m是預測值

D_1m(0)是初始值

T是體內半衰期

T2是物理半衰期(固定為175.2)

k是體內/體外放射碘的代謝比例

t是時間

• 對於每個病患，他們的D_1m(0)、T、k都不一樣，我們是否能對每個人做出一條特殊的預測線，並告訴我們什麼時候他能出院。

自訂關係式(2)

• 我們故意把Data分成兩個部分，並用第一個部分做預測式，在用此預測式來預測後面的值：

predict.func = function (D_1m.0, T1, T2 = 175.2, k, t) {
  D_1m = D_1m.0*((1-k)*exp(-log(2)*t/T1)+k*exp(-log(2)*t/T2))
  return(D_1m)
}

predict.func(D_1m.0 = 250, T1 = 3, k = 0.7, t = 0)

## [1] 250

predict.func(D_1m.0 = 250, T1 = 3, k = 0.7, t = 1)

## [1] 233.8366

predict.func(D_1m.0 = 250, T1 = 3, k = 0.7, t = 2)

## [1] 220.8678

• 我們現在再定義他的機率函數，我們假設殘差為指數常態分佈(取log後要呈現常態分佈)：

X = dat$t[1:350]
Y = dat$Y[1:350]

prob.sample = function(D_1m.0, T1, k) {
  pred.Y = predict.func(D_1m.0 = D_1m.0, T1 = T1, k = k, t = X)
  res = log(Y) - log(pred.Y)
  mean.res = 0
  sd.res = sd(res)
  log_p.res = dnorm(res, mean = mean.res, sd = sd.res, log = TRUE)
  return(-sum(log_p.res, na.rm = TRUE))
}

prob.sample(D_1m.0 = 250, T1 = 5, k = 0.1)

## [1] 575.2135

prob.sample(D_1m.0 = 200, T1 = 3, k = 0.2)

## [1] 1147.613

自訂關係式(3)

• 讓我們用最大概四估計法來取得參數吧！

fit = mle(prob.sample, start = list(D_1m.0 = 250, T1 = 10, k = 0.1), method = "SANN")
fit

## 
## Call:
## mle(minuslogl = prob.sample, start = list(D_1m.0 = 250, T1 = 10, 
##     k = 0.1), method = "SANN")
## 
## Coefficients:
##      D_1m.0          T1           k 
## 243.9800863   9.7114614   0.1224558

• 我們來比較一下預測式是否能準確預測後半段的結果(以對數後的殘差平方和為指標)：

NEW.X = dat$t[351:459]
NEW.Y = dat$Y[351:459]

pred.Y = predict.func(D_1m.0 = fit@coef[1],
                      T1 = fit@coef[2],
                      T2 = 175.2,
                      k = fit@coef[3],
                      t = NEW.X)

sum(log(NEW.Y) - log(pred.Y))^2

## [1] 124.9729

• 畫張圖看看吧！

pred.Y_1 = predict.func(D_1m.0 = fit@coef[1],
                      T1 = fit@coef[2],
                      T2 = 175.2,
                      k = fit@coef[3],
                      t = seq(0, 40, by = 0.1))

plot(dat$t, dat$Y, pch = 19, xlab = "time", ylab = "Y")
lines(seq(0, 40, by = 0.1), pred.Y_1, col = "red")

自訂關係式(4)

• 假設我們不會使用最大概似估計法，我們看了圖之後覺得時間與Y的關係是指數關係，我們試著用最簡單的方式預測：

log.Y = log(Y)

fit_log.lm = lm(log.Y~X)
pred.Y_log = fit_log.lm$coefficients[1] + NEW.X * fit_log.lm$coefficients[2]

sum(log(NEW.Y) - pred.Y_log)^2

## [1] 1.791982

• 無言了，居然比我們精心設計的關係式更準！

– 畫張圖來看看

pred.Y_2 = exp(fit_log.lm$coefficients[1] + seq(0, 40, by = 0.1) * fit_log.lm$coefficients[2])

plot(dat$t, dat$Y, pch = 19, xlab = "time", ylab = "Y")
lines(seq(0, 40, by = 0.1), pred.Y_1, col = "red")
lines(seq(0, 40, by = 0.1), pred.Y_2, col = "blue")

• 有時候以簡馭繁會不會更好?

小結

• 這節課我們學會了最大概似估計法，並且大致上已經會自訂「任何」關係式，利用最大概似估計法求解並應用來預測，但你是否注意到下列問題？

1. 關係式未必越複雜越好，事實上是越「符合事實」越好，若沒有辦法透徹了解其中關係，有時候以簡馭繁是個不錯的選擇

2. 函數「mle」在最後一個範例中，各位跑出來的結果是否都有差異呢？試著改變起始值你發現 了什麼？

• 事實上，電腦不可能會幫我們解微分，我們求得的其實是『局部極值』而非『全局極值』，但這個問題在前面的線性迴歸/邏輯斯迴歸都沒有遇到問題，這是因為他們的函數圖形較為簡單，關係式複雜時就很容易產生『局部極值』。

• 下一節課，我們將會進一步介紹在不會解微分方程時該如何求參數解，但相對應的代價就是我們通常僅能找到『局部極值』而非『全局極值』，回家想想這是重要的問題嗎?